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Abstract 

We consider the problem of interpolating an unknown multivariate 
polynomial with coefficients taken from a finite field or as numerical 
approximations of complex numbers. Building on the recent work of 
Garg and Schost, we improve on the best-known algorithm for inter- 
polation over large finite fields by presenting a Las Vegas randomized 
algorithm that uses fewer black box evaluations. Using related tech- 
^ niques, we also address numerical interpolation of sparse polynomials 

I I with complex coefficients, and provide the first provably stable algo- 

rithm (in the sense of relative error) for this problem, at the cost of 
modestly more evaluations. A key new technique is a randomization 
which makes all coefficients of the unknown polynomial distinguish- 
QQ able, producing what we call a diverse polynomial. Another departure 

\^ from most previous approaches is that our algorithms do not rely on 

root finding as a subroutine. We show how these improvements affect 
y—{ the practical performance with trial implementations. 

o 

^ 1 Introduction 

> 

Polynomial interpolation is a long-studied and important problem in com- 
puter algebra and symbolic computation. Given a way to evaluate an un- 
^ known polynomial at any chosen point, and an upper bound on the degree, 

the interpolation problem is to determine a representation for the poly- 
nomial. In sparse interpolation, we are also given an upper bound on the 
number of nonzero terms in the unknown polynomial, and the output is gen- 
erally returned in the sparse (also lacunary or supersparse) representation, 
wherein only the nonzero terms are explicitly stored. 

Applications of sparse interpolation include the manipulation and fac- 
torization of multivariate polynomials and system solving (see, e.g.. Canny 
et al. (1989); Kaltofen and Trager (1990); Diaz and Kaltofen (1995, 1998); 
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Javadi and Monagan (2007, 2009). With the advent of hybrid symbolic- 
numeric algorithms for (systems of) multivariate polynomials with approx- 
imate coefficients, we find applications of approximate sparse interpolation, 
in particular for solving non-linear systems of equations (see, e.g., Sommese 
et al. (2001, 2004); Stetter (2004)) and factoring approximate multivariate 
polynomials (see, e.g., Kaltofen et al. (2008)). 

Sparse interpolation is also a non-trivial generalization of the important 
problem of polynomial identity testing: given a black box (especially an 
algebraic circuit) computing an unknown polynomial, determine whether 
the polynomial is zero. A relevant result in our setting of sparse polynomials 
is Blaser et al. (2009); for a more in-depth discussion, we recommend the 
recent survey by Saxena (2009). 

Here we examine the sparse interpolation problem in two settings which 
have received recent attention: when the coefficients are elements of finite 
fields (particularly large finite fields, over which we have no choice) and 
when they are approximations to complex numbers. We give improvements 
over the state of the art in both of these cases, and demonstrate our new 
algorithms in practice with a full implementation in C++. 

1.1 Problem definition 

Let F be a field. A multivariate polynomial / G F[xi, . . . is said to be 
t-sparse for some t £ N if / has at most t nonzero terms in the standard 
power basis; that is, / can be written 

t 

— / ^ CjX]^ ^2 X„ 

i=l 

for coefficients Cj G F and exponent tuples {en, . . . , em ) G N" for 1 < i < t. 
If each Cj < d, then the size of this representation is 0{t) field elements plus 
0{tnlogd) bits. We seek algorithms which are polynomial-time in the size 
of this representation. 

Let / G F[3;i, . . . , x„] have degree less than d. A black box for / is 
a function which takes as input a vector (ai, . . . ,a„) G F" and produces 
/(ai, . . . , On) G F. The cost of the black box is the number of operations in 
F required to evaluate it at a given input. 

Clausen et al. (1991) showed that, if only evaluations over the ground 
field F are allowed, then for some instances at least r2(n'°^*) black box probes 
are required. Hence if we seek polynomial-time algorithms, we must extend 
the capabilities of the black box. To this end, Diaz and Kaltofen (1998) 
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introduced the idea of an extended domain black box which is capable of 
evaluating f{bi, . . . , 6„) G E for any (hi, ... , bn) G E" where E is any exten- 
sion field of F. That is, we can change every operation in the black box 
to work over an extension field, usually paying an extra cost per evaluation 
proportional to the size of the extension. 

Motivated by the case of black boxes that are division-free algebraic cir- 
cuits, we will use the following model which we believe to be fair and cover all 
previous relevant results. Here and for the remainder, M(m) is the number of 
field operations required to multiply two univariate polynomials with degrees 
less than m, and 0~{m) represents any function bounded by m(logm)*^(^). 
Using Cantor and Kaltofen (1991), M(m) G 0(mlog?nloglogm), which is 
0~{m). 

Definition 1.1. Let f G f[xi,..., x„] and i > 0. A remainder black box for 
f with size £ is a procedure which, given any monic square-free polynomial 
g G F[y] with degg = m, and any hi, . . . ,hn G F[y] with each deg/ij < m, 
produces f{hi, . . . , hn) rem g using at most I ■ M(m) operations in F. 

This definition is general enough to cover the algorithms we know of over 
finite fields, and we submit that the cost model is fair to the standard black 
box, extended domain black box, and algebraic circuit settings. The model 
makes sense over complex numbers, as we will see. 

1.2 Interpolation over finite fields 

We first summarize previously known univariate interpolation algorithms 
when F is a finite field with q elements and identify our new contributions 
here. For now, let / G Fq[a;] have degree less than d and sparsity t. We will 
assume we have a remainder black box for / with size I. Since field elements 
can be represented with 0(log q) bits, a polynomial-time algorithm will have 
cost polynomial in t, logd, and logg. 

For the dense output representation, one can use the classical method 
of Newton/Waring/Lagrange to interpolate in 0~{id) time (von zur Gathen 
and Gerhard, 2003, §10.2). 

The algorithm of Ben-Or and Tiwari (1988) for sparse polynomial in- 
terpolation, and in particular the version developed by Kaltofen and Yagati 
(1989), can be adapted to arbitrary finite fields. Unfortunately, these algo- 
rithms require t discrete logarithm computations in F*, whose cost is small if 
the field size q is chosen carefully (as in Kaltofen (2010)), but not in general. 
For arbitrary (and potentially large) q, we can take advantage of the fact 
that each discrete logarithm that needs to be computed falls in the range 
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Probes 


Probe degree 


Computation cost 


Total cost 


Dense 
Ben-Or & Tiwari 
Garg &: Schost 
Randomized G &: S 
Ours 


d 
0{t) 
0~{t^ logd) 
0~{t\ogd) 
O(logd) 


1 
1 

0~{t^ logd) 
{t^ logd) 
{t^ logd) 


0~{d) 
Oit^ + t^) 
0~{t^log^ d) 
{t^ log2 d) 
0~{t^ log2 d) 


0~{£d) 
0~{H + f + tVd) 
0~{U^ log2 d) 
0~(Vt3 log2 d) 
0~(£t2 log2 d) 



Table 1: Sparse univariate interpolation over large finite fields, 
with black box size degree d, and t nonzero terms 



[0, 1, . . . ,d - 1]. The "kangaroo method" of Pollard (1978, 2000) can, with 
high probability, compute such a discrete log with 0{^/~d) field operations. 
Using this algorithm makes brings the total worst-case cost of Ben-Or and 
Tiwari's algorithm to 0{U + f + tVd). 

The current study builds most directly on the work of Garg and Schost 
(2009), who gave the first polynomial-time algorithm for sparse interpolation 
over an arbitrary finite field. Their algorithm works roughly as follows. For 
very small primes p, use the black box to compute / modulo — 1. A 
prime p is a "good prime" if and only if all the terms of / are still distinct 
modulo — 1. If we do this for all p in the range of roughly 0{t'^ logd), 
then there will be sufficient good primes to recover the unique symmetric 
polynomial over whose roots are the exponents of nonzero terms in 
/. We then factor this polynomial to find those exponents, and correlate 
with any good prime image to determine the coefficients. The total cost is 
0~{it^ log^ d) field operations. Using randomization, it is easy to reduce this 
to 0~{£t^ log2 d). 

Observe that the coefficients of the symmetric integer polynomial in Garg 
& Schost's algorithm are bounded by 0(d*), which is much larger than the 
0(d) size of the exponents ultimately recovered. Our primary contribution 
over finite fields of size at least ^}{t'^d) is a new algorithm which avoids 
evaluating the symmetric polynomial and performing root finding over Z[?/]. 
As a result, we reduce the total number of required evaluations and develop 
a randomized algorithm with cost 0~{£t^ log^ d) , which is roughly quadratic 
in the input and output sizes. Since this can be deterministically verified in 
the same time, our algorithm (as well as the randomized version of Garg & 
Schost) is of the Las Vegas type. 

The relevant previous results mentioned above are summarized in Ta- 
ble 1, where we assume in all cases that the field size q is "large enough". 
In the table, the "probe degree" refers to the degree of g in each evaluation 
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of the remainder black box as defined above. 
1.3 Multivariate interpolation 

Any of the univariate algorithms above can be used to generate a multivari- 
ate polynomial interpolation algorithm in at least two different ways. For 
what follows, write p{d, t) for the number of remainder black box evaluations 
required by some univariate interpolation algorithm, A{d,t) for the degree 
of the remainder in each evaluation, and ip{d, t) for the number of other field 
operations required besides black box calls. Observe that these correspond 
to the first three columns in Table 1. 

The first way to adapt a univariate interpolation algorithm to a mul- 
tivariate one is Kronecker substitution: given a remainder black box for 
an unknown / G ^[xi, . . . , Xn ], with each partial degree less than d, we 
can easily construct a remainder black box for the univariate polynomial 
/ = /(x, x'^, x'^^ , . . . , x'^" ^) G F[x], whose terms correspond one-to-one with 
terms of /. This is the approach taken for instance in Kaltofen (2010, §2) 
for the interpolation of multivariate polynomials with rational coefficients. 
The cost is simply the cost of the chosen underlying univariate algorithm, 
with the degree increased to d". 

The other method for constructing a multivariate interpolation algo- 
rithm is due to Zippel (1990). The technique is inherently probabilistic and 
works variable-by-variable, at each step solving a number of txt transposed 
Vandermonde systems, for some t < t. Specifically, each system is of the 
form Ax = b, where A is a t x t matrix of scalars from the coefficient field 
¥q. The vector v consists of the output of t remainder black box evalua- 
tions, and so its elements are in Fg[y], and the system must be solved modulo 
some g £ as specified by the underlying univariate algorithm. Observe 

however that since A does not contain polynomials, computing x = A~^b 
requires no modular polynomial arithmetic. In fact, using the same tech- 
niques as Kaltofen and Yagati (1989, §5), employing fast dense bivariate 
polynomial arithmetic, each system can be solved using 



field operations. 

Each transposed Vandermonde system gives the remainder black box 
evaluation of each of t univariate polynomials that we are interpolating in 
that step. The number of such systems that must be solved is therefore 
p{d,t), as determined by the underlying univariate algorithm. Finally, each 
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Kronecker 


Zippel 


Dense 




O'iintd) 


Ben-Or & Tiwari 


O'iit + t2 + t(f^/2) 


0~{nt^ + nt^Vd + int^) 


Garg &: Schost 


0~{£n'^t^ log2 d) 


0~{£nt^ log2 d) 


Randomized G & S 


{inH^ log2 d) 


{£nt'^ log2 d) 


Ours 


O {£nH^ log2 d) 


0~(£nt3 log2 d) 



Table 2: Sparse multivariate interpolation over large finite fields, 
with black box size n variables, degree d, and t nonzero terms 

of the t univariate interpolations proceeds with the given evaluations. The 
total cost, over all iterations, is 

0~{lnt- /\{d,t)- p{d,t)) 

field operations for the remainder black box evaluations, plus 

0~{nt^{d, t) + Int ■ A{d, t)) 

field operations for additional computation. Zippel (1990) used the dense 
algorithm for univariate interpolation; using Ben-Or and Tiwari's algorithm 
instead was studied by Kaltofen and Lee (2003). 

Table 2 summarizes the cost of the univariate algorithms mentioned 
above applied to sparse multivariate interpolation over a sufficiently large 
finite field, using Kronecker's and Zippel's methods. 

For completeness, we mention a few more results on closely related prob- 
lems that do not have a direct bearing on the current study. Grigoriev et al. 
(1990) give a parallel algorithm with small depth but which is not com- 
petitive in our model due to the large number of processors required. A 
practical parallel version of Ben-Or and Tiwari's algorithm has been devel- 
oped by Javadi and Monagan (2010). Kaltofen ct al. (1990) and Avendaho 
et al. (2006) present modular algorithms for interpolating polynomials with 
rational and integer coefficients. However, their methods do not seem to 
apply to finite fields. 

1.4 Approximate Polynomial Interpolation 

In Section 4 we consider the case of approximate sparse interpolation. Our 
goal is to provide both a numerically more robust practical algorithm, but 
also the first algorithm which is provably numerically stable, with no heuris- 
tics or conjectures. We define an "e-approximate black box" as one which 
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evaluates an unknown t-sparse target polynomial / G C[x], of degree d, with 
relative error at most e > 0. Our goal is to build a t-sparse polynomial g 
such that II/ — g'||2<e||/||2- A bound on the degree and sparsity of the tar- 
get polynomial, as well as e, must also be provided. In Section 4 we formally 
define the above problem, and demonstrate that the problem of sparse in- 
terpolation is well-posed. We then adapt our variant of the Garg and Schost 
(2009) algorithm for the approximate case, prove it is numerically accurate 
in terms of the relative error of the output, and analyze its cost. We also 
present a full implementation in Section 5 and validating experiments. 

Recently, a number of numerically-focussed sparse interpolation algo- 
rithms have been presented. The algorithm of Gicsbrccht et al. (2009) is 
a numerical adaptation of Ben-Or and Tiwari (1988), which samples / at 
0{t) randomly chosen roots of unity a; G C on the unit circle. In particular, 
u is chosen to have (high) order at least the degree, and a randomization 
scheme is used to avoid clustering of nodes which will cause dramatic ill- 
conditioning. A relatively weak theoretical bound is proven there on the 
randomized conditioning scheme, though experimental and heuristic evi- 
dence suggests it is much better in practice. Cuyt and Lee (2008) adapt 
Rutishauser's qd algorithm to alleviate the need for bounds on the partial 
degrees and the sparsity, but still evaluate at high-order roots of unity. Ap- 
proximate sparse rational function interpolation is considered by Kaltofen 
and Yang (2007) and Kaltofen et al. (2007), using the Structured Total 
Least Norm (STLN) method and, in the latter, randomization to improve 
conditioning. Approximate sparse interpolation is also considered for inte- 
ger polynomials by Mansour (1995), where a polynomial-time algorithm is 
presented in quite a different model from ours. In particular the evaluation 
error is absolute (not relative) and the complexity is sensitive to the bit 
length of the integer coefficients. 

Note that all these works evaluate the polynomial only on the unit circle. 
This is necessary because we allow and expect / to have very large degree, 
which would cause a catastrophic loss of precision at data points of non-unit 
magnitude. Similarly, we assume that the complex argument of evaluation 
points is exactly specified, which is again necessary because any error in the 
argument would be exponentially magnified by the degree. 

The primary contribution of the work in this paper is to provide an 
algorithm with both rigorously provable relative error and good practical 
performance. Our algorithm typically requires 0~{t^ log^ d) evaluations at 
primitive roots of unity of order 0~{t^ log d) (as opposed to order d in pre- 
vious approaches). We guarantee that it finds a t-sparse polynomial g such 
that 11^ — / II2 < 2e||/||2. An experimental demonstration of the numerical 
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robustness is given in Section 5. 



2 Sparse interpolation for generic fields 

Here and for the remainder, we say a polynomial / is t-sparse if it can be 
written as a sum of at most t nonzero coefficients times a monomial. We 
assume the unknown polynomial / is always univariate. This is without 
loss of generality, as we can use the Kronecker substitution as discussed 
above. The exponential increase in the univariate degree only corresponds 
to a factor of n increase in logdeg /, and since our algorithms will ultimately 
have cost polynomial in logdeg/, polynomial time is preserved. 

Assume a fixed, unknown, t-sparse univariate polynomial / G f[x] with 
degree at most d. We will use a remainder black box for / to evaluate 
/rem(xP — 1) for small primes p. We say p is a "good prime" if the sparsity 
of frem{xP — 1) is the same as that of / itself — that is, none of the 
exponents are equivalent modulo p. 

The following lemma shows the size of primes required to randomly 
choose good primes with high probability. 

Lemma 2.1. Let f G f[x] be a t-sparse polynomial with degree d, and 
let A = max (21, Wt{t — 1) Ind]). A prime chosen at random in the range 
A, . . . , 2A is a good prime for f with probability at least 1/2. 

Proof. Write ei, . . . ,et for the exponents of nonzero terms in /. If p is a bad 
prime, then p divides {ej — et) for some i < j. Each ej — ei < d, so there 
can be at most log^ d = In d/ In A primes that divide each ej — Cj. There are 
exactly (2) such pairs of exponents, so the total number of bad primes is at 
most l)ln(i)/(21nA). 

From Rosser and Schoenfeld (1962, Corollary 3 to Theorem 2), the total 
number of primes in the range A, . . . , 2A is at least 3A/(51n A) when A > 
21, which is at least t{t — l)lnd/lnA, at least twice the number of bad 
primes. □ 

Now observe an easy case for the sparse interpolation problem. If a 
polynomial / G F[x], has all coefficients distinct; that is, / = X]i<j<tC«^^* 
and Ci = Cj ^ i = j, then we say / is diverse. To interpolate a diverse 
polynomial / G ^[x], we first follow the method of Garg and Schost (2009) 
by computing /rem(xP* — 1) for "good primes" pi such that the sparsity of 
/rem(xP' — 1) is the same as that of /. Since / is diverse, /rem(3;P' — 1) is 
also diverse and in fact each modular image has the same set of coefficients. 



8 



Algorithm 1: Generic interpolation 



Input: /i G M>o, T,D,q G N, and a remainder black box for unknown 

T-sparse / E F[x] with deg/ < D 
Output: t e N, ei, . . . , G N, and ci, . . . , q € F such that 
/ — J2i<i<t Cia;*^* 

1 t^O 

2 A ^ max(21, [|r(r- l)lnD]) 

3 for [log2(3//u)] primes p e {A, . . . , 2A} do 
Use black box to compute fp = f{x) rem (re'' — 1) 
if fp has more than t terms then 

6 t -^r- sparsity of fp 

7 Q ^ V 

s a element of F such that Pr[/(aa;) is not diverse] < /u/3 
9 gg f{ax) rem(a;^ — 1) 

10 ci , . . . , Q ■(— nonzero coefficients of 

11 ei,...,et^O 

12 for [21n(3//i) +4(lnL>)/(lnA)] primes p e {X, ... ,2X} do 

13 Use black box to compute gp = f{ax) iem{x^ — 1) 

14 if gp has exactly t nonzero terms then 

15 for i = 1, . . . ,t do Update with exponent of Cj in gp modulo 
p via Chinese remaindering 



16 for i = 1, . . . ,t do Cj ^ CiO 

17 return f{x) = Ei<i<tCia;''' 



Using this fact, we avoid the need to construct and subsequently factor the 
symmetric polynomial in the exponents. Instead, we correlate like terms 
based on the (unique) coefficients in each modular image, then use simple 
Chinese remaindering to construct each exponent from its image modulo 
each Pi. This requires only 0(log d) remainder black box evaluations at good 
primes, gaining a factor of t improvement over the randomized version of 
Garg and Schost (2009) for diverse polynomials. 

In the following sections, we will show how to choose an a G F so that 
f{ax) — which we can easily construct a remainder black box for — is 
diverse. With such a procedure. Algorithm 1 gives a Monte Carlo algorithm 
for interpolation over a general field. 
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Theorem 2.2. With inputs as specified, Algorithm 1 correctly computes the 
unknown polynomial f with probability at least 1 — fj,. The total cost in field 
operations (except for step 8) is 

Proof. The for loop on hne 3 searches for the true sparsity t and a single 
good prime g. Since each prime p in the given range is good with probability 
at least 1/2 by Lemma 2.1, the probability of failure at this stage is at most 
/i/3. 

The for loop on line 12 searches for and uses sufficiently many good 
primes to recover the exponents of /. The product of all the good primes 
must be at least D, and since each prime is at least A, at least (lnL')/(lnA) 
good primes are required. 

Let n = [21n(3/^) + 4(lnD)/(lnA)] be the number of primes sampled in 
this loop, and k = [(lnZ))/(lnA)] the number of good primes required. We 
can derive that (n/2 — A;)^ > (ln(3//i) + fc)^ > (n/2) ln(3//u), and therefore 
exp(— 2(| — ky/n) < jj/S. Using Hoeffding's Inequality (Hoeffding, 1963), 
this means the probability of encountering fewer than k good primes is less 
than /i/3. 

Therefore the total probability of failure is at most fi. For the cost 
analysis, the dominating cost will be the modular black box evaluations in 
the last for loop. The number of evaluations in this loop is 0(log(l//i) + 
(log Z))/ (log A)), and each evaluation has cost 0{£ • M(A)). Since the size of 
each prime is 0((log£')/(logT + loglog£')), the complexity bound is correct 
as stated. □ 

In case the bound T on the number of nonzero terms is very bad, we 
could choose a smaller value of A based on the true sparsity t before line 8, 
improving the cost of the remainder of the algorithm. 

In addition, as our bound on possible number of "bad primes" seems to 
be quite loose, a more efficient approach in practice would be to replace the 
for loop on line 12 with one that starts with a prime much smaller than A 
and incrementally searches for the next larger primes until the product of 
all good primes is at least D. We could choose the lower bound to start 
searching from based on lower bounds on the birthday problem. That is, 
assuming (falsely) that the exponents are randomly distributed modulo p, 
start with the least p that will have no exponents collide modulo p with high 
probability. This would yield an algorithm more sensitive to the true bound 
on bad primes, but unfortunately gives a worse formal cost analysis. 



10 



3 Sparse interpolation over finite fields 



We now examine the case that the ground field F is the finite field with q 
elements, which we denote ¥q. First we show how to effectively diversify 
the unknown polynomial / in order to complete Algorithm 1 for the case of 
large finite fields. Then we show how to extend this to a Las Vegas algorithm 
with the same complexity. 



3.1 Diversification 

For an unknown / G ¥q[x] given by a remainder black box, we must find 
an a so that f{ax) is diverse. A surprisingly simple trick works: evaluating 
f{ax) for a random nonzero a £¥q. 

Theorem 3.1. Forq > T(T—1)T+1 and any T- sparse polynomial f G ^glx] 
with deg/ < D, if a is chosen uniformly at random from ¥*, the probability 
that /(ax) is diverse is at least 1/2. 

Proof. Let t < T be the exact number of nonzero terms in /, and write 
/ = l^i<j<j Cja;^% with nonzero coefficients Cj G F* and ei < 62 < • • • < e^. 
So the ith coefficient of f{ax) is Qa*^'. 

If f{ax) is not diverse, then we must have c^a'^' = Cja^^ for some i 7^ j. 
Therefore consider the polynomial A G ¥q[y] defined by 

We see that f{ax) is diverse if and only if A(a) 7^ 0, hence the number of 
roots of A over Fg is exactly the number of unlucky choices for a. 

The polynomial A is the product of exactly (2) binomials, each of which 
has degree less than D. Therefore 

^ ^ TiT-l)D 
deg A < ^ , 

and this also gives an upper bound on the number of roots of A. Hence 
q — 1 > 2deg^, and at least half of the elements of F* are not roots of A, 
yielding the stated result. □ 

Using this result, given a black box for / and the exact sparsity t of f, 
we can find an a G Fg such that f{ax) is diverse by sampling random values 
a G Fq, evaluating f{ax) rem — 1 for a single good prime p, and checking 
whether the polynomial is diverse. With probability at least 1 — n, this 
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will succeed in finding a diversifying a after at most [log2(l/^)] iterations. 
Therefore we can use this approach in Algorithm 1 with no effect on the 
asymptotic complexity. 



3.2 Verification 

So far, Algorithm 1 over a finite field is probabilistic of the Monte Carlo 
type; that is, it may give the wrong answer with some controllably-small 
probability. To provide a more robust Las Vegas probabilistic algorithm, we 
require only a fast way to check that a candidate answer is in fact correct. 
To do this, observe that given a modular black box for an unknown T- 
sparse / G ^q[x] and an explicit T-sparse polynomial g G Fg[2;], we can 
construct a modular black box for the 2T-sparse polynomial f — g of their 
difference. Verifying that f = g thus reduces to the well-studied problem of 
deterministic polynomial identity testing. 

The following algorithm is due to Blaser et al. (2009) and provides this 
check in essentially the same time as the interpolation algorithm; we restate 
it in Algorithm 2 for completeness and to use our notation. 



Algorithm 2: Verification over finite fields 



Input: T, D,q £ N and remainder black box for unknown T-sparse 

/ G ¥g[x] with degf <D 
Output: ZERO iff / is identically zero 

1 for the least (T — 1) log2 D primes p do 

2 Use black box to compute fp = /rem(x^ — 1) 

3 if /p / then return NONZERO 

4 return ZERO 



Theorem 3.2. Algorithm 2 works correctly as stated and uses at most 

O [IT log • M (T log D • (log T + loglog D))) 
field operations. 

Proof. For correctness, notice that the requirements for a "good prime" 
for identity testing are much weaker than for interpolation. Here, we only 
require that a single nonzero term not collide with any other nonzero term. 
That is, every bad prime p will divide ej — ei for some 2 < j < T. There can 
be log2 D distinct prime divisors of each Cj — ei , and there are T — 1 such 
differences. Therefore testing that the polynomial is zero modulo — 1 for 
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the first (T — 1) log2 D primes is sufficient to guarantee at least one nonzero 
evaluation of a nonzero T-sparse polynomial. 

For the cost analysis, the prime number theorem (Bach and Shallit, 1996, 
Theorem 8.8.4), tells us that the first (T — 1) log2 D primes are each bounded 
by 0(r -logD- (logT + loglogD)). The stated bound follows directly. □ 

This provides all that we need to prove the main result of this section: 

Theorem 3.3. Given q > T(T — 1)D + 1, any T,D £N, and a modular black 
box for unknown T-sparse f £ ^q[x] with degf < D, there is an algorithm 
that always produces the correct polynomial f and with high probability uses 
only 0~ (^IT"^ log^ Z?) field operations. 

Proof. Use Algorithms 1 and 2 with = 1/2, looping as necessary until the 
verification step succeeds. With high probability, only a constant number 
of iterations will be necessary, and so the cost is as stated. □ 

For the small field case, when q E 0{T^D), the obvious approach would 
be to work in an extension E of size 0(log T + log D) over Fg. Unfortunately, 
this would presumably increase the cost of each evaluation by a factor of 
logD, potentially dominating our factor of T savings compared to the ran- 
domized version of Garg and Schost (2009) when the unknown polynomial 
has very few terms and extremely high degree. 

In practice, it seems that a much smaller extension than this is sufficient 
in any case to make each gcd(ej — ei,q — 1) small compared to g — 1, but we 
do not yet know how to prove any tighter bound in the worst case. 

4 Approximate sparse interpolation algorithms 

In this section we consider the problem of interpolating an approximate 
sparse polynomial / G C[x] from evaluations on the unit circle. We will 
generally assume that / is t-sparse: 



We require a notion of size for such polynomials, and define the coefficient 



/= E 



CjX^' for Ci £ C and ei < ■ ■ ■ < et = d. 



(4.1) 



i<j<t 



2-norm of / = ^| 



0<i<d 



fi X BjS 
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The following identity relates the norm of evaluations on the unit circle 
and the norm of the coefficients. As in Section 2, for / G C[x] is as in (4.1), 
we say that a prime p is a good prime for / if p f (cj — ej) for all i ^ j. 

Lemma 4.1. Let f G C[x], p a good prime for f, and a; G C a pth primitive 
root of unity. Then 

^ 0<i<p 

See Giesbrecht and Roche (2010, Theorem 2.9). 

We can now formally define the approximate sparse univariate interpo- 
lation problem. 

Definition 4.2. Let e > and assume there exists an unknown t-sparse 
f G C[x] of degree at most D. An e-approximate black box for f takes an 
input ^ G C and produces a 7 G C such that I7 — f{(,)\ < e|/(OI- 

That is, the relative error of any evaluation is at most e. As noted 
in the introduction, we will specify our input points exactly, at (relatively 
low order) roots of unity. The approximate sparse univariate interpolation 
problem is then as follows: given D, T G N and 5 > e > 0, and an e- 
approximate black box for an unknown T-sparse polynomial / G C[x] of 
degree at most D, find a T-sparse polynomial g £ C[x] such that ||/ — (7II2 ^ 

^Ilfi'll2- 

The following theorem shows that t-sparse polynomials are well-defined 
by good evaluations on the unit circle. 

Theorem 4.3. Let e > and f G C[x] be a t-sparse polynomial. Suppose 
there exists a t-sparse polynomial g £ C[x] such that for a prime p which is 
good for both f and f — g, and pth primitive root of unity uj £ C, we have 

1/(0;*) - 5(^^)1 <e|/(a;*)| for < i < p. 

Then \\f — g\\2 ^ ^\\f\\2- Moreover, if go G C[x] is formed from g by deleting 
all the terms not in the support of f, then ||/ — (^o|l2 ^ 2e II/II2. 

Proof. Summing over powers of uj we have 

5^ \f{co^)-g{u:'^)\'<e' ^ \f{co^)\\ 

0<i<p 0<i<p 

Thus, since p is a good prime for both f — g and /, using Lemma 4.1, 
P • 11/ -5II2 < - P- II/II2 and 11/ - C/II2 < e ||/||2- 
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Since g — go lias no support in common with /, 

II5-50II2 < II/-5II2 < e|l/ll2- 

Thus 

II/-50II2 = 11/ -5 + (5-50)112 

< 11/ -5II2 + 115 -50II2 < 2e II/II2 • □ 

□ 

In other words, any t-sparse polynomial whose values are very close to / 
must have the same support except possibly for some terms with very small 
coefficients. 

4.1 Computing the norm of an approximate sparse polyno- 
mial 

Let < e < 1/2 and / £ C[x] a i-sparse polynomial for which we are given 
an e-approximate black box. We first consider the problem of computing 

II/II2- 



Algorithm 3: Approximate norm 
Input: T,D gN and e-approximate black box for unknown T-sparse 

/ G C[x] with degf <D 
Output: o" G M, an approximation to II/II2 

1 A ^ max(21, l)lnd]) 

2 Choose a prime p randomly from {A, . . . , 2A} 

3 uj exp{2-iri/p) 

i w •'r- {/{io^), . . . , f{uj^^^)) e computed using the black box 
5 return (l/y/p) ■ \\w\\2 



Theorem 4.4. Algorithm 3 works as stated. On any invocation, with prob- 
ability at least 1/2, it returns a value a G M>o such that 

(l-2e)||/||2 <a< (l + e)||/||2. 

Proof. Let v = {f{u^), . . . , f{ioP~^)) G be the vector of exact evaluations 
of /. Then by the properties of our e-approximate black box we have w = 
V + eA, where |Aj| < |/(a;*)| ioi < i < p, and hence || AII2 < ||f ||2- By the 
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triangle inequality \\w\\2 < ||^^||2 + e || A||2 < (l + e)||f||2. By Lemmas 2.1 
and 4.1, ||u||2 = -v/p||/|l2 with probability at least 1/2, so (1/y^) • \\w\\2 < 
(1 + e) II/II2 with this same probability. 

To establish a lower bound on the output, note that we can make er- 
ror in the evaluation relative to the output magnitude: because e < 1/2, 
|/(a;*) — Wi\ < 2e\wi\ for < i < p. We can write v = w + 2eV, where 
IIVII2 < ||w||2- Then ||v||2 < (1 + 2e) \\w\\2, and (1 - 2e) II/II2 < {l/^/p) ■ 
\\w\\2. □ 

4.2 Constructing an e-approximate remainder black box 

Assume that we have chosen a good prime p for a t-sparse / G F[x]. Our goal 
in this subsection is a simple algorithm and numerical analysis to accurately 
compute / rem — 1. 

Assume that / rem x^ — 1 = Ylo<i<p exactly. For a primitive pth root 
of unity u G C, let V{uj) £ C^^p be the Vandermonde matrix built from the 
points l,w, . . .,ujP-^. Recahthat V{uj)-{bo,.. • ,6p-i)^ = • • . 

and V{uj~^) = p-V{ijj)~^ . Matrix vector product by such Vandermonde ma- 
trices is computed very quickly and in a numerically stable manner by the 
Fast Fourier Transform (FFT). 



Algorithm 4: Approximate Remainder 
Input: An e-approximate black box for the unknown t-sparse 

f £ C[x], and p G N, a good prime for / 
Output: h G C[x] such that || (/rem x^ — 1) — h\\2 < e ||/||2- 

1 w ^ {f{uj^), . . . , /{ujP^^)) G C^' computed using the 

e-approximate black box for / 

2 n (1/p) ■ V{uj~^)w G C'P using the FFT algorithm 

3 return h = Y^Q^i^pUiX^ G C[x] 



Theorem 4.5. Algorithm 4 works as stated, and 

II (/rem j;^ — 1) — /i||2<e||/||2. 

It requires 0{plogp) floating point operations and p evaluations of the black 
box. 

Proof. Because / and / rem — 1 have exactly the same coefficients {p is 
a good prime for /), they have exactly the same norm. The FFT in Step 
2 is accomplished in 0{plogp) floating point operations. This algorithm 
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is numerically stable since • V{uj^^) is unitary. That is, assume 

V = {f{ujQ),...,f{ojP''^)) G CP is the vector of exact evaluations of /, so 
11^^ — w^|l2 — ^11^112^7 the black box specification. Then, using the fact that 



\/p||/ll2' 

II (/ rem x^~^] 



1 



1 



V{uj-'^) -{v-w) 



-V{u-^)v--V{lo-^)w 
p p 

1 „ 



V — w 



< 



□ 



4.3 Creating e-diversity 

First, we extend the notion of polynomial diversity to the approximate case. 

Definition 4.6. Let f G C[x] be a t-sparse polynomial as in (4.1) and 
6 > € > such that |cj| > 5 II/II2 for 1 < i <t. The polynomial f is said to 
be e-diverse if and only if every pair of distinct coefficients is at least e II/II2 
apart. That is, for every 1 < i < j < t, \ci — Cj\ > e \\f\\2- 

Intuitively, if (e/2) corresponds to the machine precision, this means that 
an algorithm can reliably distinguish the coefficients of a e-diverse polyno- 
mial. We now show how to choose a random a to guarantee e-diversity. 

Theorem 4.7. Let (5 > e > and f £ C[x] a t-sparse polynomial whose non- 
zero coefficients are of magnitude at least S \\f\\2- If s is a prime satisfying 
s > 12 and 

t{t-l)<s< 3.1-, 
e 

then for = e^'^'/* an s-PRU and G N chosen uniformly at random from 
{0, 1, . . . , s — 1}, fiC^x) is e-diverse with probability at least ^. 

Proof. For each 1 < i < t, write the coefficient q in polar notation to base 
C as Cj = ri(^% where each rj and 6i are nonnegative real numbers and 
n > S\\f\\2- 

Suppose f{C,^x) is not e-diverse. Then there exist indices 1 < i < j < t 
such that 



Because minfrj,rj] 



> 5\\f\ 



5II/II2- C 



hfcei 
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the value of the left hand side is at least 



. Dividing out Q 



-kci 



, we get 



^k(e 



-ei) 



- 5 
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By way of contradiction, assume there exist distinct choices of k that 
satisfy the above inequahty, say ki,k2 G {0, . . . , s — 1}. Since and 
(^^j-^i are a fixed powers of C not depending on the choice of k, this means 



Because s is prime, Cj ^ ej, and we assumed ki 7^/02, the left hand side is 
at least |^ — 1|. Observe that 2tt/s, the distance on the unit circle from 1 
to C, is a good approximation for this Euclidean distance when s is large. 
In particular, since s > 12, 



and therefore |C — 1| > 6^/2{V^ — l)/s > 6.2/s, which from the statement 
of the theorem is at least 2e/6. This is a contradiction, and therefore the 
assumption was false; namely, there is at most one choice of k such that the 
i'th and j'th coefficients collide. 



We note that the diversification which maps /(x) to f{C^x) and back is 
numerically stable since ^ is on the unit circle. 

In practice, the previous theorem will be far too pessimistic. We there- 
fore propose the method of Algorithm 5 to adaptively choose s, 6, and 
simultaneously, given a good prime p. 

Suppose there exists a threshold 5 G N such that for all primes s > S, 
a random sth primitive root of unity makes f{C^x) e-diverse with high 
probability. Then Algorithm 5 will return a root of unity whose order is 
within a constant factor of S, with high probability. Prom the previous 
theorem, if such an S exists it must be 0{t'^), and hence the number of 
iterations required is O(logt). 

Otherwise, if no such S exists, then we cannot diversify the polynomial. 
Roughly speaking, this corresponds to the situation that / has too many 
coefficients with absolute value close to the machine precision. In this case, 
we can simply use the algorithm of Garg and Schost (2009) numerically, 
achieving the same stability but using a greater number of evaluations and 
bit operations. It is possible to establish an adaptive hybrid between our 
algorithm and that of Garg and Schost (2009) by making / as e-diverse as 
possible given our precision. The non-zero coefficients of / are clustered into 
groups which are not e-diverse (i.e., are within e||/||2 of each other). We 




6 



C-l| ^ ^/2(^/3-l) /2 
27r/s 7r/6 
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Algorithm 5: Adaptive diversification 





Input: e-approximate black box for /, known good prime p, known 






sparsity t 




Output: C, k such that f{C^x) is e-diverse, or FAIL 


1 


s 


I, 6 ^ oo, fp<—0 


2 


while s < t'^ and i^{coeffs c of fs s.t. \c\ > 5} < t do 


3 




s ^ least prime > 2s 


4 




C ^ exp(27ri/s) 


5 




k ^ random integer in {0, 1, . . . , s — 1} 


6 




Compute fs = fiC'^x) rem — 1 


7 




6 ^ least number s.t. all coefficients of fs at least S in absolute 






value are pairwise e-distinct 


8 


if 5 > 2e then return FAIL 


9 


else return 



can use the symmetric polynomial reconstruction of Garg and Schost (2009) 
to extract the exponents within each group. 

4.4 Approximate interpolation algorithm 

We now plug our e-approximate remainder black box, and method for mak- 
ing / e-diverse, into our generic Algorithm 1 to complete our algorithm for 
approximate interpolation. 

Theorem 4.8. Let 5 > 0, / G C[x] with degree at most D and sparsity 
at most T, and suppose all nonzero coefficients c of f satisfy \ c\ > 5\\f\\2- 
Suppose also that e < 1.56/{T{T — 1)), and we are given an e-approximate 
black box for f. Then, for any fi < 1/2 we have an algorithm to produce a 
g G C[x] satisfying the conditions of Theorem 4-3. The algorithm succeeds 
with probability at least 1 — /U and uses 0~(T^ ■ log(l//i) • log^ D) black box 
evaluations and floating point operations. 

Proof. Construct an approximate remainder black box for / using Algo- 
rithm 4. Then run Algorithm 1 using this black box as input. On step 8 
of Algorithm 1, run Algorithm 5, iterating steps 5~7 [log2(3/^)] times on 
each iteration through the while loop to choose a diversifying q = with 
probability at least 1 — fi/3. 

The cost comes from Theorems 2.2 and 4.5 along with the previous 
discussion and Theorem 4.7. □ 



19 



Observe that the resulting algorithm is Monte Carlo, but could be made 
Las Vegas by combining the finite fields zero testing algorithm discussed in 
Section 3.2 with the guarantees of Theorem 4.3. 

5 Implementation results 

We implemented our algorithms in C++ using the GNU Multiple Precision 
Arithmetic Library (GMP, http : //gmplib . org/) and Victor Shoup's Num- 
ber Theory Library (NTL, http://www.shoup.net/ntl/) for the exponent 
arithmetic. For comparison with the algorithm of Garg and Schost (2009), 
we also used NTL's squarefree polynomial factorization routines. We note 
that, in our experiments, the cost of integer polynomial factorization (for 
Garg & Schost) and Chinese remaindering were always negligible. 

In our timing results, "Determ" refers to the deterministic algorithm as 
stated in Garg and Schost (2009) and "Alg 1" is the algorithm we have 
presented here over finite fields, without the verification step. We also de- 
veloped and implemented a more adaptive, Monte Carlo version of these 
algorithms, as briefly described at the end of Section 2. The basic idea is to 
sample modulo — 1 for just one prime p G 0(t^log(i) that is good with 
high probability, then to search for much smaller good primes. This good 
prime search starts at a lower bound of order 0(f^) based on the birthday 
problem, and finds consecutively larger primes until enough primes have 
been found to recover the symmetric polynomial in the exponents (for Garg 

6 Schost) or just the exponents (for our method). The corresponding im- 
proved algorithms are referred to as "G&S MC" and "Alg 1++" below, 
respectively. 

Table 3 summarizes some timings for these four algorithms over the finite 
field Z/65521Z. This modulus was chosen for convenience of implementa- 
tion, although other methods such as the Ben-Or and Tiwari algorithm 
might be more efficient in this particlar field since discrete logarithms could 
be computed quickly. However, observe that our algorithms (and those from 
Garg and Schost) have only poly-logarithmic dependence on the field size, 
and so will eventually dominate. 

The timings are given in seconds of CPU time on a 64-bit AMD Phenom 
II 3.2GHz processor with 512K/2M/6M cache, compiled using GCC 4.4.3 
with the -03 flag. Note that the numbers listed reflect the base-2 logarithm 
of the degree bound and the sparsity bound for the randomly-generated test 
cases. 

The timings are mostly as expected based on our complexity estimates. 
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0.31 


0.01 


0.01 
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3.66 


1.06 


0.65 


20 
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9.32 
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Table 3: Finite Fields Algorithm Timings 
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Noise 


Mean Error 


Median Error 


Max Error 





4.440 e-16 


4.402 e-16 


8.003 e-16 


±10-12 


1.113e-14 


1.119e-14 


1.179 e-14 


±10-9 


1.149 e-11 


1.191 e-11 


1.248 e-11 


±10-^ 


1.145 e-8 


1.149 e-8 


1.281 e-8 



Table 4: Approximate Algorithm Stability 

and also confirm our suspicion that primes of size 0{t^) are sufficient to 
avoid exponent collisions. It is satisfying but not particularly surprising to 
see that our "Alg 1±±" is the fastest on all inputs, as all the algorithms 
have a similar basic structure. Had we compared to the Ben-Or and Tiwari 
or Zippel's method, they would probably be more efficient for small sizes, 
but would be easily beaten for large degree and arbitrary finite fields as their 
costs are super-polynomial. 

The implementation of the approximate algorithm uses machine double 
precision (IEEE), the built-in C±± complex<double> type, and the popu- 
lar Fastest Fourier Transform in the West (FFTW, http : //www . f f tw . org/) 
package for computing FFTs. Our stability results are summarized in Ta- 
ble 4. Each test case was randomly generated with degree at most 2^^ and 
at most 50 nonzero terms. We varied the precision as specified in the table 
and ran 10 tests in each range. Observe that the error in our results was 
often less than the e error on the evaluations themselves. 

Both implementations are released under an MIT-style licence and are 
available from the second author's website at 
http : //www . cs . uwaterloo . ca/~droche/ diverse/. 

6 Conclusions 

We have shown how to use the idea of diversification to improve the complex- 
ity of sparse interpolation over large finite fields by a factor of i, the number 
of nonzero terms. We achieve a similar complexity for approximate sparse 
interpolation, and provide the first provably numerically stable algorithm 
for this purpose. Our experiments confirm these theoretical results. 

Numerous open problems remain. A primary shortcoming of our algo- 
rithms is the quadratic dependence on t, as opposed to linear in the case of 
dense interpolation or even sparse interpolation in smaller or chosen finite 
fields using the Ben-Or and Tiwari algorithm. It seems that reducing this 
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quadratic dependency will not be possible without a different approach, be- 
cause of the birthday problem embedded in the diversification step. In the 
approximate case, a provably numerically stable algorithm for sparse inter- 
polation with only 0{t) probes is still an open question. And, while general 
backward error stability is not possible in the high degree case, it would be 
interesting in the case of low degree and many variables. 
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